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Abstract 

We performed a high statistics simulation of Ising spins coupled to 
2D quantum gravity on toroidal geometries. The tori were triangulated 
using the Regge calculus approach and contained up to 512^ vertices. 
We used a constant area ensemble with an added interaction term, 
employing the dl/l measure. We find clear evidence that the critical 
exponents of the Ising phase transition are consistent with the static 
critical exponents and do not depend on the coupling strength of the 

interaction term. We definitively can exclude for this type of mo- 
del a behaviour as predicted by Boulatov and Kazakov [Phys. Lett. 
B186, 379 (1987)] for Ising spins coupled to dynamically triangulated 
surfaces. 
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1 Introduction 



With the advances made in string theory there was a rising interest in 2D 
quantum gravity, motivated by the fact, that a string moving in N dimension 
is equivalent to N scalar fields coupled to 2D quantum gravity [1]. For N 
equal to zero we then deal with pure gravity. An important step forward 
was made when Kazakov [2] suggested a model of Ising spins living on the 
vertices of (f>^ graphs that was solvable, and in its dual form, equivalent to 
Ising spins coupled to a dynamical triangulated surface (DTS). It turned out 
[3] that the set of critical exponents of the Ising transition is very different 
from the static Onsager critical exponents (see Table I below). Some time 
later Knizhnik et al. (KPZ) [4] found the same set of critical exponents with 
methods of conformal field theory for matter of central charge c = 1/2, which 
is supposed to be the continuum limit of Kazakov's model. The exponents 
have been confirmed in a variety of numerical studies [5-8] on cfr" graphs as 
well as with the dual method of dynamical triangulated surfaces. 

One of the oldest methods to study general relativity numerically was 
suggested by Regge [9]. He uses a simplicial approximation to a manifold, 
where the underlying lattice has fixed coordination number, and takes the 
link lengths as gravitational degrees of freedom, whereas in the DTS method 
it is just the opposite, the coordination number varies and the simplices have 
fixed link lengths. The relation of Regge calculus to classical relativity has 
been extensively studied [10, II], and its continuum limit is reasonably well 
understood. There has been up to now only one simulation [12] with rather 
low statistics of Ising spins coupled to gravity using the Regge approach, 
which suggests, that the critical exponents are the same as in the fiat case. 
This is somewhat disturbing, because it is generally believed, that both ap- 
proaches should describe the same model in the continuum limit. We found 
it therefore worthwhile to investigate the Regge approach again on larger lat- 
tices and with a better statistics, to confirm or disprove the results obtained 
in Ref. [12]. 

2 The model and simulation techniques 

We simulated the gravitational interaction using the Regge calculus where 
the underlying manifold is discretized with a fixed triangulation and the link 
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lengths are subject to variations. We used the usual transcription [13, 14] of 
continuum quantities like the metric g and the scalar curvature R into the 
Regge approach, namely 



'd'x^^) (1) 

i 

j (Px^l^)R{x) 2J2S,. (2) 

i 

jSx^)R'{x) ^4^^, (3) 



where 8i is the deficit angle at the vertex i defined as 

,5, = 27r- (4) 

all t sharing i 

and 9i{t) is the dihedral angle associated with the triangle t. The area Ai is 
taken to be the baricentric area associated with the site z, 

A = E^^t' (5) 

where At denotes the area of the triangle t, which is one of several popu- 
lar choices that are believed to be equivalent in the continuum limit. We 
simulated the partition function 

Z = J2 I Dii{l) exp (-/(/) - KE{1, s)) (6) 

M 

where /(/) is defined as the gravitational action 

/(/) = ^(AA. + a|), (7) 

i ' 

and 

E{hs) = \ E M"-^r (8) 

edges lij 

is the energy of Ising spins s^, Si = ±1, which are located at the vertices i of 
the lattice. Here the volume Aij associated with a link lij is defined as 

^ E l^^- (9) 

triangles t D l^j 
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In two dimensions the Einstein-Hilbert action (2) is in virtue of the Gauss- 
Bonnet theorem proportional to the Euler characteristic x, and therefore a 
topological invariant which does not contribute to the gravitational action 
(7). In higher than two dimensions a interaction term is sometimes added 
to guarantee the boundedness of the gravitational action from below. This 
is not necessary in two dimensions where we included such a term to probe 
its influence on the continuum limit. 

The form of the path integral measure Dfi[l) in (6) is already not clear 
in the continuum formalism. The most simple and most often used choice 
on the lattice is Dfi[l) = JJdl/l^ which we also adopted here. We simulated 
the gravitational action using the standard single-hit Metropolis update. In 
addition to the usual Metropolis procedure a change in link length is only 
accepted, if the links of a triangle fulfill the triangle inequality. As Ising 
update we used the single-cluster (Wolff) algorithm [15] which prevents the 
critical slowing down near the phase transition. Between measurements we 
performed n = 2 ... 4 Monte Carlo steps consisting of one lattice sweep to 
update the link lengths lij followed by a single-cluster fiip to update (a frac- 
tion of) the spins Si. We checked in some cases that varying the relative 
frequency of link and spin updates does not change the results within error 
bars. 

We simulated the partition function (6) on triangulated tori of size N = 
with fixed coordination number q = 6. This gives rise to 2N triangles and 
3A^ link variables. The principal simulations were performed at a = 0.001 and 
the couplings K = 1 and K = 1.025 for L = 6,8,10,12,16,32,64,100,128, 
200, 256, and 512. Additional simulations were performed with a = and 
O.I at K = 1.025, using lattices of size L = 8,16,32,64,100,128 and 256. 
Because of the scale invariance of the measure we could rescale each link 
when proposing a link update such that the total area was kept fixed to its 
initial value A = N. The difference of the model defined by (6) and the 
Ising model on a static triangular lattice is that the spins are coupled by 
geometric weight factors Wij = Ai^ jP- which can fiuctuate around the static 
value = a/3/6, with a peak at Wij < w^^ and a long tail for large Wij. 
With decreasing coupling a the peak sharpens and its location shifts to small 
values of Wiy 

For each run we recorded the time series of the energy density e = E j A, 
the magnetization density m = AiSij A and the Liouville field density 
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= Xlj liiAi/A. For each lattice we performed about 50000 measurements. 
From an analysis of the time series we found integrated autocorrelation times 
for the energy and the magnetization of about 1—4 (in units of measurements) 
for all lattice sizes. To obtain results for the various observables O at K values 
in an interval around the simulation point Kq^ we applied the reweighting 
method [16]. Since we recorded the time series this amounts to computing 

with AK = K — Kq. To obtain errors we devided each run into 20 blocks 
and computed standard Jackknife errors. At a = 0.001 where we had two 
simulations at different K values, we combined the results according to their 
errors [17, 18]. 

From the time series we computed the Binder parameter [19], 

It is well known that the Ul{K) curves for different L cross around (A^c, U*) 
with slopes oc L^^'^^ apart from confluent corrections explaining small syste- 
matic deviations. This allows an almost unbiased estimate of the critical 
coupling A^c, the critical correlation length exponent and the renormalized 
charge U* . The slopes can be conveniently calculated as 

dUL , A , , (m^E) im^E)\ , , 

— ^ = l-[/i A -2Y^ + V^ • 12 
dK I y^ ) y^ ) J 

We further analyzed the (finite lattice) susceptibility, 

x(A0 = A((m2)-(|m|)2), (13) 
the susceptibility in the disordered phase, 

X'(A0 = A((m2)), (14) 

the specific heat, 

C{K) = K'A{{e') - {en (15) 
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and studied the (finite lattice) magnetization at its point of inflection, (|m|) |inf . 
The inflection point can be obtained from the maximum of d(\m\) j dK , which 
can be calculated as 

'^-Ml = ^E){\m\)-{E\m\). (f6) 

Further useful quantities are the logarithmic derivatives 

d\n{\m\) {E\m\) 

and 



Another gravitational quantity of interest is the Liouville field ^{x) = In ^ 
In the discretized version its lattice average reads as = l/A^lilnAi, and 
the associated lattice Liouville susceptibility is defined as 

X^L) = A{{^') - i^f). (19) 



3 Results 

By applying reweighting techniques we first determined the maxima of x, 
C, d{\m\) / dK ^ dln{\m\) / dK ^ and dln{m^) / dK . The location of the ma- 
xima provided us with five sequences of pseudo-transition points Kina.x{L) 
for which the scaling variable x = {K^a^^{L) — Kc)L^^'^ should be constant. 
Using this information we then have several possibilities to extract the criti- 
cal exponent p from (linear) least square fits of the finite-size scaling (FSS) 
Ansatz^ dUL/dK = L^^^foix) or dln{\m\P)/dK = L^^^fpix) to the data at 
the various K^a_x{L). The fit range was chosen such that the goodness-of-fit 
parameter Q was always above 0.1. 

For the very extensive simulations at a = 0.001 all values are in good 
agreement with the Onsager value p = 1 at a 2 % level, giving rise to an 

"'^By writing the scaling Ansatz as = N^l^" instead of = L^l" we trivially get the same 
estimates for 'IjvD. This would be directly comparable with DTS simulations where D is 
possibly a non-trivial internal fractal dimension; cp. Table 1. 
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average of Ijv = 1.00(1). For the other two couphngs, a = 0.1 and a = 0, 
the data scatter a bit more but their averages Ijv = 0.98(l)(a = 0.1) and 
Ijv = 0.95(2)(a = 0) are still compatible with v = 1. 

Assuming thus v = 1 we have next determined estimates for Kc from 
the Binder parameter crossings and the scahng of the various K^a_x{L). The 
crossings of the curves Ul{K) with L and L' approach Kc as 

A-x = + - 1), (20) 

where b = L'/L^ k, is a constant, and confluent corrections are neglected. 
This method, valid for large 6, turned out to be the most precise one, leading 
for a = 0.001 to an estimate of the critical coupling Kc = 1.02652(6). For 
a = 0.1 we obtain by the same method Kc = 1.0292(1), and for a = we 
find Kc = 1.0230(2). Performing linear least square fits with p = 1 to the 
various pseudo-transition points K^a_x{L) gave consistent estimates of Kc for 
all five sequences. As final value we quote their average Kc = 1.02650(8), 
where the error is taken as the minimum among the five estimates. This is 
probably an overestimate, but since the five sequences are of course correlated 
a more refined optimization of the final estimate would require a careful 
analysis of cross-correlations which we have not performed here. For the other 
couplings a = 0.1 and a = 0.0 we obtain Kc = 1.0301(2), and Kc = 1.0243(4), 
respectively, again in reasonable agreement with the estimate from the Binder 
parameter crossings. Combining this information we use in further analyses 

Kc = 1.0234 ± 0.0002 (a = 0.0), (21) 
Kc = 1.0265 ± 0.0001 (a = 0.001), (22) 
Kc = 1.0295 ± 0.0001 (a = 0.1). (23) 

In particular we can now extract p also from the scaling of dU / dK and 
dln{\m\''') / dK at Kc^ in good agreement with our previous results. 

Another quantity of interest is the asymptotic limit U* of the Binder 
parameter at Kc- We first looked at the scaling of = U{K^) for fixed 
b = L' I L = 2, where one can assume a scaling of the form 

[/^(A^x) = [/* + civ-'". (24) 

From a three-parameter fit we obtain U* = 0.616(9) and u = 0.6(5) for 
a = 0.001. Another somewhat more stable way to obtain an estimate for U* 
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Figure 1: Finite-size scaling of the Binder parameter values Ul{K) for a = 
0.001 and K = 1.0265 ~ Kc- The solid line is a three-parameter fit Ul{K) = 
U* + cL-"^, yielding U* = 0.6117(14) and uj = 1.1(4). 

is by fitting the values of Ul at Kc according to (24). For a = 0.001 we obtain 
from the three-parameter fit shown in Fig. 1 U* = 0.6117(14) and io = 1.1(4). 
The uncertainty in A^c, however, implies a further error of AU* ~ 0.0035, 
so that we finally quote U* = 0.612(5). From three-parameter fits at Kc 
we obtain for a = 0.1 U* = 0.615(6) and uj = 0.7(4). For a = 0.0 the 
errors turned out to be larger, leading to an estimate of U* = 0.609(33) and 
UJ = 0.3(7). 

These values are practicably indistinguishable from estimates for the regu- 
lar square lattice, U* = 0.615(10) [20] and U* = 0.611(1) [21], or Poissonian 
random lattices, U* = 0.615(7) [22]. This gives further evidence that the cri- 
tical behaviour of Ising spins on Regge lattices is governed by the standard 
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(Onsager) universality class. We are not aware of any estimates for U* in the 
DTS approach. 

To extract the critical exponent ratio 7/;/ we used the scaling x — 
L^'^'^fsi^) the previously discussed points of constant x, as well as the 
scaling of x' at Kc- The quality of the fits for Xmax can be inspected in Fig. 2. 
As final results we find -f/p = 1.735(5)(a = 0.001), -f/p = 1.728(3)(a = 0.1), 
and 7/;/ = 1.717(6)(a = 0). The values for 7/;/ for the different values of a 
are compatible with each other, but are all slightly below the Onsager value 
of 7/;/ = 1.75. Due to their respective error range, however, they are still 
consistent with the fiat space exponent ratio. 

To extract the magnetical critical exponent ratio jS / p we used that (|m|) = 
L~^^'^ fi{x) at all constant x-values. Another method is too look at the scaling 
of d{\m\)ldK = L^^-f^y^fsix). Because the errors on the different estimates 
turned out to vary over a large range we chose to compute error-weighted 
averages. Using our average values for p we obtain the final estimates of 
/3/p = 0.127(3) (a = 0.001), /3/p = 0.123(2) (a = 0.1), and /3/p = 0.123(4) 
(a = 0.0). Again we see little infiuence of the curvature square term, and the 
results are again in agreement with the Onsager result fi j p = 0.125. 

Having found so far overwhelming evidence for the Onsager universality 
class, we expect also the specific-heat exponent a to take on the Onsager 
value, namely a = 0. In this case we expect a logarithmic divergence like 

C{x,L) = A{x) + B{x)\nL. (25) 

Indeed the data at the different fixed values of x can all be fitted nicely with 
this Ansatz. In particular, for a fit of Cmax with 11 data points at a = 0.001 
we obtain A = 0.17(1), B = 0.369(5), with Q = 0.80. Of course, we also did 
an unbiased 3-parameter fit of the form 

C(x, L) = A{x) + B{x)L"l\ (26) 

which gave us in the case of the fit of Cmax with a = 0.001 and 8 data points 
A = 7.9(3.8), B = -7.9(3.7), and ajp = -0.06(5), with a quality Q = 0.12, 
compare Fig. 3. For comparison, we included also the best fit with Ansatz 
(25), and a linear least-square fit with the DTS prediction ajp = —2/3. 
For the other two a-values we get A = 7.4(8.1), B = —7.4(7.8), and ajp = 
-0.07(9), with a quality Q = 0.08 for a = 0.1, and A = 7(12), B = -7(12), 
and ajp = —0.06(13), with a quality Q = 0.32 for a = 0.0. Similar fits of C at 
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1 2 3 4 5 6 7 



In L 



Figure 2: Double logarithmic finite-size scaling plot of the susceptibility ma- 
xima Xmax for a = 0.0,0.001, and 0.1. To disentangle the curves we added 
an offset of —2 (2) to the data for a = 0.0 (a = 0.1). The slopes are in all 
three cases compatible with the Onsager value 7/;/ = 1.75 for regular static 
lattices. 
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the other Kina.x{L) sequences as well as at Kc gave consistent results. Another 
method, which in previous studies [23] gave improved estimates over (26), is 
to fit the energy density according to e[L) = a-\- BL^"~^'^^'^. Three-parameter 
fits at gave {a - I) j v = -0.98(5)(a = 0.001), (a - l)lv = -0.98(4)(a = 
0.1), and (a — 1) j v = — 1.06(4)(a = 0.0). Overall we can conclude that also 
for aj V our data is fully consistent with the static Onsager value of zero. 

We found in all FSS analyses that the added interaction term did not 
affect the FSS behaviour. To get a final estimate we therefore computed 
a weighted average of the three simulations with different coupling a. The 
results, put into comparison with the KPZ predictions, can be found in Ta- 
ble 1. We used the usual scaling relations rj = 2 — ^ j v and (5 = l-|-7//3, as 
well as our measured value for v to convert the exponent ratios. 

Let us finally consider the critical behaviour of the Liouville field which 
is the interesting variable in the gravitational sector. The FSS prediction is 
Xif{L) = iv^^~''*'-'/6(x) [12]. Using the 7 largest values for L in our simulation 
for a = 0.001 at Kc a linear-least-square fit yields (2 — rj,^) = 0.047(28), with 
Q = 0.15, and for a = 0.1 we get with 6 data points (2 — rj,^) = 0.040(13), 
with Q = 0.79. For a = 0, on the other hand, we find with 7 data points 
(2 — T],^) = 2.01(6), with Q = 0.19. The observed pronounced dependence of 
T],^ on a is plausible because a positive a term tends to suppress curvature 
fiuctuations and stabilizes in this way also the area fiuctuations. Therefore 
one should only consider the case with a = 0, which gives an observed value 
of T],^ ^ that is consistent with a massless free field behaviour of (p. 

4 Concluding remarks 

We have performed a fairly detailed high statistic study of the Ising model 
coupled to quantum gravity via the Regge calculus approach. Using the path 
integral measure JJdl/l we have found that an included interaction term 
showed no effect on the nature of the phase transition, and that the critical 
exponents of the Ising transition still agree with the Onsager exponents for 
regular static lattices. If one believes in KPZ scaling then Regge calculus 
can survive as a tool to probe quantum gravity only, if one can modify it 
to reproduce the KPZ results. The two most likely alterations would be 
to implement a different coupling of the Ising spins to gravity, or to use a 
different, maybe more physically motivated, measure. In fact, exploratory 
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Figure 3: Finite-size scaling of the specific-heat maxima for a = 0.001. Also 
shown are a logarithmic fit Cmax = A -\- BlnL^ a power-law fit Cmax = 
A+BL"^\ as well constrained power-law fit assuming the DTS prediction 
a/v = -2/3. 
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simulations [24] with a different measure seem to imply that the choice of 
measure plays a more dominant role than has been thought so far. 
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a 


/3 


7 




1] 


V 


DTS 


-1 


0.5 


2 


5 


2/3* 


1.5* 


Onsager 





0.125 


1.75 


15 


0.25 


1 


this MC study 


^ 


0.126(2) 


1.75(2) 


14.9(3) 


0.272(3) 


1.01(1) 



Table 1: Comparison of our Monte Carlo results with the exact results for 
the Ising model on static lattices (Onsager) and the results of Boulatov and 
Kazakov [3] for the Ising model on dynamical graphs (DTS). The values 
marked with a star were computed from hyperscaling relations with D = 2, 
thereby neglecting possible scaling effects due to the internal fractal dimen- 
sion in the DTS approach. 
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